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Abstract 



O 

A hybrid Monte Carlo (HMC) approach is employed to quantify the influence of 

inelastic deformation on the microstructural evolution of polycrystalline materials. 

This approach couples a time explicit material point method (MPM) for deformation 

'a 

with a calibrated Monte Carlo model for grain boundary motion. A rate-independent 
*t3 ■ crystal plasticity model is implemented to account for localized plastic deformations 

in polycrystals. The dislocation energy difference between grains provides an addi- 



tional driving force for texture evolution. This plastic driving force is then brought 



J> ■ into a MC paradigm via parametric links between MC and sharp-interface (SI) ki- 

OV 
O ■ netic models. The MC algorithm is implemented in a parallelized setting using a 

in 

checkerboard updating scheme. As expected, plastic loading favors texture evolution 

l> ■ 

O ■ for grains which have a bigger Schmid factor with respect to the loading direction, 

o: 

and these are the grains most easily removed by grain boundary motion. A macro- 
scopic equation is developed to predict such texture evolution. 

Keywords : Monte Carlo, grain boundary, plasticity, anisotropy, texture, driving 
force 



Corresponding author. Tel.: 718-260-3082; fax: 718-260-3532. E-mail: lzhang@poly.edu (L. 
Zhang) 

1 



1. Introduction 



The prediction of microstructural evolution in response to thermo-mechanical 
loading is important for materials design, processing or thermomechanical fatigue 
phenomena. Computational modeling of evolving texture in response to large plas- 
tic deformation and recrystallization has been studied extensively |l|, y, y|, Ifl but 
less so than that produced by thermally-induced stresses-i.e. stress-induced texture 
evolution [5|, |6j. We consider a thermo-mechanical setting in which temperature 
changes cause stresses to develop due to geometrical constraints. The temperature 
is sufficiently high to generate grain boundary motion and yet low enough such that 
recrystallization does not occur. The induced stresses may be associated with both 
elastic and plastic deformation [5fl. 

n 

In a previous work, a Hybrid Monte Carlo (HMC) approach [7[ was developed by 
combining a MC algorithm for grain boundary motion |8| with the Material Point 
Method (MPM) [9J for elastic deformation. Purely elastic driving forces, originating 
from the anisotropic mechanical response of individual grains, are treated as a bulk 



body force in a Potts model for grain boundary evolution [lfj. The approach is time 



accurate through the use of parametric links 



to sharp-interface (SI) kinetics 11] 



It also takes advantage of the 



of the driving force 
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act that MC grain boundary mobility is independent 



uence of inelastic de- 
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texture evolution 



The present work extends this paradigm to include the infli 
formation on texture evolution [6|. As in the elastic study 

is assumed to be dominated by grain boundary kinetics 15|, [16|. Furthermore, we 
consider infinitesimal deformation to distinguish the stress-induced texture from de- 
formation texture. The latter is associated with grain and lattice rotation in response 



to finite deformations 13, ll8[. A stochastic, crystal plasticity model, developed from 
rate- independent crystal plasticity 19], is applied within the MPM framework as 
the constitutive model to capture the elasto-plastic response of a poly crystalline me- 
dia 20|. As opposed to conventional deterministic algorithms [l9[, the stochastic 
algorithm relies on a MC routine to determine the activated slip system which is 
therefore referred to as the Monte Carlo Plasticity (MCP). When plastic deforma- 
tion occurs, dislocations are generated, stored and annihilated within the microstruc- 
ture. The heterogeneous distribution of these dislocations within the polycrystalline 
medium constitutes a plastic driving force for grain boundary migration. This is 
treated as a body force within the MC kinetics using parametric links between MC 
and SI models. A Red/Black (RB) updating scheme is used to parallelize the MC 
algorithm 21], although other methods might also be useful 22J. This parallelized 
HMC approach is used to investigate the microstructural evolution of nickel poly- 
crystals under plastic loading. As expected, the grains with smaller Schmid factors 
gradually dominate the polycrystalline system. The data is subsequently used to 
construct a macroscopic kinetic equation to predict the evolution of microstructure. 

2. Monte Carlo crystal plasticity 

Plastic response of polycrystalline materials is treated thr oug h a classical rate- 
independent small deformation crystal plasticity formulation 23]. The foundations 
of the constitutive model assume that the elasto-plastic response of single crystals is 



dominated by slip deformation mechanisms [24], |25|, |26j. A successful numerical al- 
gorithm must carry out three tasks: the determination of activated slip systems; the 
calculation of the plastic slip on each activated slip system; and, the solution of redun- 

methods 
have been devised and successfully implemented in deterministic formats 19l. |27I|. 



dant constraints associated with a hardening law [191.127]. Various numerica 



As opposed to deterministic algorithms, the current work adopts a probabilistic 
approach borrowed from concepts in statistical mechanics in which only one slip 
system is activated during each time step. Plastic slip is therefore treated as a series of 
discrete, probabilistic events that mimic the sequential accumulation of dislocations 
at the lattice scale. This Monte Carlo crystal plasticity (MCP) is algorithmically 
simple because plastic slip can be resolved through the solution of one equation with 
no redundant constraints. On the other hand, the associated computational steps 
has to be sufficiently small such that a sequence of single slips mimics multiple slip 
behavior. A probabilistic algorithm, detailed in what follows, is used to determine 
which slip system is chosen at each step. The constitutive framework and stress 
updating routine are otherwise standard 19]. 

Given a set of potentially activated slip systems £ = (1,2, ...,ri), identified through 
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the elastic energy of a 



comparison of resolved shear stress with slip resistance 

crystal, E a , can be calculated if each slip system of the set C is individually activated. 

This generates n possible states for the deformed crystal. The probability, p a , of a 



slip system being selected is computed using the partition function [29], Z : 

n 
Z = V e -0E« 

Pa 



a=l 
„-0E n 



:i) 



where a is the index of a potentially activated slip system, and /3 is the inverse of 
the fundamental temperature. Dislocation energy can be ignored in Eqn. flTJ due to 
the fact that an isotropic hardening model is used 19j, 
the set £ is activated when the following criterion is met: 



28] . The ith slip system of 



Y,Pa~Pi< R< Y,Pa- 



(2) 



«=1 



«=1 



Here R is taken as a random number between and 1. With the activated slip 
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Elastic moduli 
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123.5 GPa 
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8 MPa 


Saturation stress 


Ts 


0.122 GPa 


Initial hardening 


ho 


0.44 GPa 


Inverse fundamental temp 


P 


1000.0 



Table 1: Nickel properties required in MCP model [l9l[28J. 

system determined, the deformation can easily be parsed into elastic and plastic 
components. 

As a verification of this MCP method, a uniaxial tension test was carried out on 
a single crystal of nickel. The material properties 
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28| required for this model 



are listed in Tabled] A prescribed strain increment, Ae x , was enforced at each time 
step: 



Ae 3 



e 











-¥ 











-¥ 



(3) 



with e = 1.0 x 10~ 6 . Fig. [T]compares the results obtained using the MCP logic with 
their counterparts obtained with a deterministic model wherein a singular value 
decomposition (SVD) algorithm is applied to solve the ill-conditioned constraint 
equations. The two approaches clearly agree. 

3. Grain boundary kinetics 



Sharp-interface, grain boundary kinetics ll| is the basis for our grain boundary 
kinetics model. It can be developed by considering a bi-crystal subjected to uniaxial 
loading that results in plastic deformation as shown in Fig. [2j The dislocation 
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Figure 1: Comparison of stress-strain relationship between MCP and SVD algorithms with a pre- 
scribed strain increment. The discrete dots and solid line represent MCP and SVD results, respec- 
tively. 
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Figure 2: A single grain boundary with SI notation. 



energy, E p , and dislocation density, p are related using the relation [30] 



E P = \ P Gb% P=(^) 2 - (4) 

Here G is the shear modulus, b g is the Burgers vector, and o$ is the flow stress. The 



thermodynamic driving force for interfacial accretion is given by 11] 



/, = \E p \ + k 7 * , (5) 

where 7* is the capillary driving force, k is the local mean curvature, and [~] is the 
jump in a field across the grain boundary. The elastic driving force is temporarily 
suppressed in order to study the influence of the plastic behavior on texture evolu- 
tion. The Herring relation 3l| is used to describe the accretive normal speed of the 
interface, v n : 

v n = m s f s , (6) 

where m s is grain boundary mobility. 

The MC paradigm is intended to implement SI kinetics within a computationally 
efficient setting. Continuum, deterministic fields are linked to discrete, probabilistic 
counterparts, and the physical domain is discretized into a square lattice. A Q-state 
Potts model [lOj is used to represent a crystalline system and grain orientation is 
described by an integer- valued spin field, qt, which ranges from 1 to Q. The system 
Hamiltonian is then written as 

N M N 

j=l n=l i=l 

where bi is the spatially varying bulk energy associated with the i th lattice, J quqn is 
the interaction energy between the neighboring spin fields qi and q n , M is the number 
of neighbors considered on the selected lattice, and S is the Kronecker delta. Kinetics, 

7 



within the MC paradigm, is translated to a series of probabilistic trials carried out 
at all lattice sites. The acceptance of a trial event is determined by a probability 
function based on the change in the Hamiltonian associated with a trial flip. Plastic 
effects can therefore be accounted for if the plastic driving force is treated as a bulk 
energy term in Eqn. (J7j). 

Parametric links between SI and MC models are used to convert SI parameters 
into MC format [7J: 

- 7*^0 

Is 



A 2 Jl 2 Q 
2bE 

-Ml 



m. 



A 4 m(cw,#K 



J T mc Eoto 

where 7% and m s are reasonable grai n boundary stiffness and grain boundary mo- 



bility obtained from the literature 
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33], while b s is the computed bulk energy 



density associated with the SI model. The counterparts in the MC paradigm, which 
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361 ] , b and m 
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can be analytically or numerically derived, are 7 
respectively. A characteristic energy, E , length, l , and time, to along with the MC 
lattice size A, interaction energy J, effective temperature a mc , inclination angle 6 
and time step r mc , are also required to bridge the SI and MC paradigms [13]. These 
links endow the MC simulation with physical time and length scales. 



4. A Parallelized Monte Carlo Algorithm 

In a typical MC implementation, an update of variables immediately follows the 
acceptance of a trial event as in the Metropolis algorithm for example 8|, 137]. To 
improve the efficiency of this approach within a parallelized setting, we adopted 
a Red/Black (RB) updating rule |2l| wherein the domain is decomposed into a 



checkerboard, and each MC step is divided into red/black and black/red half steps 



3N 



In contrast to the 



This is an alternative to the N-fold way algorithm [22 . 
large number of updates associated with a standard MC time step, states are not 
updated in the RB format until each half step finishes. When applied within a parallel 
environment, the domain is divided uniformly among processors with synchronous 
communication between processors at each half-step. 



The approach is first implemented within a 2D, bi-crystal setting |l3j. A circular 
grain is placed at the center of a square domain and the inner grain shrinks or 
grows due to the combined effects of capillary and bulk driving forces. The two- 
state square lattice Ising model is then employed to represent the physical system. 
MC simulations were carried out on a 500 x 500 grid, at an effective temperature 
a mc = 1.4. The evolution of inner grain was simulated with four combinations 
of capillary and bulk driving forces. Capillary driving force was fixed by giving 
a constant interaction energy, i.e., J = 1, while bulk driving force, which is half 
of the bulk energy difference per site between inner and outer grains, varies from 
—0.1 to 0.1. In particular, the bulk energy term will be used later on to represent 
mechanical effect (dislocation energy) in realistic plastic deformation. Corresponding 
SI simulations were also run using the established parametric links between MC and 
SI models. Fig. [3] (a) indicates good agreement between the two models in the 
prediction of inner grain area as a function of time. An additional test case was also 
performed for a polycrystalline system where the Q-state (Q = 91) Potts model with 
isotropic interaction energy (J = 1) was used. Many MC simulations were carried 
out at effective temperatures a mc = 1.0 and a mc = 0.5, respectively. Fig. |3j(b) shows 
the measured average grain size with linear fits to the data. The linear grain growth 



behavior is consistent with that of the Isotropic Grain Growth Theory 40|, |41|. In 



later applications, a relatively high MC effective temperature is preferred to remove 
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Figure 3: Validation of the parallelized MC algorithm in 2D: (a) Comparison of the inner grain area 
evolution of a bi-crystal between MC and SI results with four combinations of capillary and bulk 
driving forces at a mc ~lA. Solid line and discrete points represent MC and SI results respectively. 
Grid size is 500 x 500. (b) Isotropic grain growth at two MC temperatures for 2-D polycrystals. 
The gray and black dots demonstrate a rnc =0.5 and a mc =1.0, respectively. Grid size is 500 x 500. 
MC results are averaged over 10 samples. 



lattice pinning 40] 



5. Texture evolution 

The aforementioned parallelized MC algorithm was next applied within the pre- 
vious polycrystalline setting to study the influence of bulk energy distribution on 
the evolution of texture. The interaction energy was taken to be isotropic and a 
Gaussian distribution for the bulk energy was assigned to each orientation. Fig. [4] 
(a) shows the orientation-dependent Gaussian shape bulk energy: 



b(n) 



-Exp[- 



(n - 46)' 



(9) 



(320°- 5 ri 2 x 46 2 J ' 
where n is orientation, and b is the bulk energy associated with this material. The 
bulk energy is adimensioned within a pure MC setting. MC simulations were carried 
out at an effective temperature a mc = 1.0. Fig. |4] (b) presents the orientation 
distribution at three time slices. As physically expected, materials points with a 
lower bulk energy are favored. This indicates that grains which experience the least 
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Figure 4: The impact of anisotropic bulk energy on polycrystalline texture evolution: (a) Gaussian 
distribution of orientation-dependent nondimcnsional bulk energy, and (b) orientation distributions 
at three time slices. The grid size is 500 x 500, and a mc = 1.0. Simulation results were averaged 
over 100 runs. 



plastic deformation will be favored as opposed to the ones have larger dislocation 
content. 

The parallelized HMC approach was then used to consider the microstructural 
evolution of three-dimensional, polycrystalline nickel in response to plastic deforma- 
tion. Attention was restricted to texture development in response to plastic deforma- 
tion which is dominated by grain boundary kinetics rather than that of triple junction 



kinetics 



15 



)l.ll6l|. To meet this requirement, all numerical experiments were performed 
at an MC effective temperature (a mc = 2.6), which is a non-physical temperature. 
Note that the corresponding physical temperature (as opposed to MC temperature) 
is still not high enough to initiate recrystallization. In addition, elastic driving forces 
were suppressed for the sake of clarity. The MCP framework is implemented in MPM 
as a plastic constitutive model mediating the mechanical response 19[ . The required 
fitting parameters are listed in Table [TJ Material properties at room temperature 
were adopted even though they may be temperature dependent 42[. In order to il- 
lustrate the methodology and facilitate the interpretation of the results, grains were 
distinguished by a single angle of rotation with respect to the z-direction and allowed 
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it to vary from 0° to 90° in 1° angle steps. The other two Euler angles were held 
fixed. A similar setting has been previously consider with deformation restricted to 
the elastic regime [7|]. 

As mentioned previously, plastic slip occurs when the resolved shear stress on 
a slip system exceeds its threshold resistance. As a consequence, orientations with 
bigger Schmid factor will be plastically deformed first. Fig. [5] shows the orientation 
dependent Schmid factor of a face centered cubic (fee) single crystal. Orientations 
near 22.5° and 67.5° have the largest Schmid factor. These orientations slip more 
easily and thus accumulate dislocations the fastest. Therefore, subsequent grain 
boundary motion will tend to remove such grains. In the current study, an 20 MPa 
uniaxial loading was applied in the x-direction in 100 loading increments. At each 
step, the strain increment associated with each MPM particle is computed using 
the MPM algorithm, and the MCP algorithm is then called to partition elastic and 
plastic deformations. The quasi-static state of the polycrystalline system is reached 
after a series of iterations when either the relative stress increment or relative strain 
increment meets the prescribed convergence criterion for all the particles [7J. After 
each MCP step, dislocation energy is computed and converted into the MC domain 
using our parametric links. The microstructure is subsequently evolved for one MC 
step, equal to 10 minutes J7|. Figs. El(a) and (b) describe the evolution of orientation 
distribution at five different time slices. Negative rotation angles are used to represent 
orientation from the 46° to 90°. As expected, texture evolution favors grains which 
have smaller Schmid factors with respect to the loading axis. 

Interestingly, orientations between 23° and 45° are not selected with the same 
frequency as orientations between 0° and 22° even though the two groups have iden- 
tical Schmid factors. This is because the effective Young's moduli of the former 
group elements are larger than that of their counterparts. Therefore less elastic, and 
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Figure 5: Orientation-dependent Schmid factors. Each orientation is described by a single orienta- 
tion angle with respect to the z-direction. 
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Figure 6: Texture evolution under a uniaxial plastic loading in the x-direction. Each orientation 
is represented by a single Euler angle. Simulation data are shown as discrete points, which are 
averaged over 10 independent runs. 



thus more plastic strain, results from the same loading. To illustrate this princi- 
ple more clearly, pure mechanical simulations were performed in single crystals of 
various orientations. Fig. [7J (a) shows that more plastic deformation is obtained 
on the 22° orientation than the other two orientations. Fig. [7J (b) highlights the 
differences in the mechanical response between 0° and 45° orientations which have 
the same Schmid factor. As expected, more plastic deformation is observed in the 
45° orientation. 
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Figure 7: Anisotropic behaviors of single crystals under plastic deformation: (a) stress-strain rela- 
tionship for three single crystals with different Schmid factors, (b) stress-strain relationship for two 
crystals with the same Schmid factors. 



(U = 



0.06 
0.05 
0.04 
0.03 
0.02 
0.01 









x 






18 

CD 

u 
£ 16 

jg 14 
g 12 

i io 

to 

U 8 


















/ i 


* \ ■■ 


1000 Mins 
800 Mins 
600 Mins 


















u? 




>\\ 
























•*5a 








41e -0.0185t _ 






\*~ 








o^*3t 
















O^ 



































































-20 -10 10 20 

Rotation Angle (Deg) 



400 



600 800 

Time (mins) 



1000 



Figure 8: Numerical fits of texture evolution with plastic influences: (a) Simulation data are shown 
as discrete points along with Gaussian fits as solid curves, and (b) Evolution of the Gaussian 
variance for the orientation distributions. 



To quantify the effect of texture evolution with plasticity, the texture histograms 
were fitted to a time dependent Gaussian equation |7J]: 

1 A2 



f(M 



Exp[- 



(10) 



g(t)y/2ir iL 2g(t)^ 
where / is the orientation probability density, <p is the orientation angle, and g is a 
time- dependent Gaussian variance. The fitted results are shown as solid curves in 
Fig. [8] (a) along with the Gaussian variance described in Fig. (b). 
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6. Discussion and Summary 

The effects of plastic deformation on texture evolution are quantified with a pre- 
viously developed hybrid algorithm that blends a discrete, probabilistic algorithm 
for grain boundary motion with a continuum, deterministic model of elastic defor- 
mation. Plastic deformation is accounted for using a Monte Carlo plasticity model 
in which slip occurs through a sequence of probabilistically determined, single slip 
events. Mechanical loading results in an inhomogeneous distribution of dislocation 
energy which amounts to a plastic driving force for texture evolution. Grains with 
less damage grow at the expense of those which have been more heavily deformed. 

In the current approach, the dislocation density associated with the area swept 
by a moving interface is replaced with that of the overtaking grain. This implies 
that, when a grain boundary moves, the dislocation content adjacent to the grain 
boundary is extended into the material through which the grain boundary migrates. 
However, this is not always the case in reality. On the other hand, if one assume that 
no dislocations were transmitted with the boundary motion, then a new dislocation 



free grain would be introduced and the material would recrystallize [30|. Though 
recrystallization undoubtedly occurs, and is certainly important in the materials 
studied and at these levels of deformation, it is not the only option. A grain boundary 
can simply migrate in a heavily deformed polycrystal |43|. The migration of a grain 
boundary in a deformed polycrystal involves the propagation of some, but not always 
all, of the dislocation content that are adjacent to it. The authors are not aware of 
a meaningful quantitative description of how the dislocation structure propagates 
with a migrating grain boundary, and thus assume that nearly all the dislocations 
are carried with a boundary when it moves. 

This hybrid computational methodology enables the temporal evolution of the 
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microstructure under thermomechanical loading. It offers an alternative to sharp- 
interface and phase-field modeling. The resulting texture evolution maps are ex- 
pected to be useful in materials and process design, thermomechanical fatigue as 
well as in part performance assessment throughout service life. The current work 
was limited to a single angle of misorientation between grains. An extension of this 
work to include fully arbitrary misorientation is underway, which will allow us to cap- 
ture the influence of plastic deformation on the texture evolution of more realistic 
material systems. 
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